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A VARIABLE-STEP  CENTRAL  DIFFERENCE  METHOD  FOR  STRUCTURAL  DYNAMICS  ANALYSIS 

PARTI:  THEORETICAL  ASPECTS 


K.  C.  PARK  AND  P.  G.  UNDERWOOD 

Applied  Mechanic*  Laboratory 
Lockheed  Palo  Alto  Reaearch  Laboratory 
Palo  Alto,  California  94104 


ABSTRACT 

Three  major  theoretical  aspecta  of  variable-step 
explicit  integration  procedures  are  proposed  and  ana- 
lyzed, These  Include  basic  fixed-step  integration 
formulas  for  no  damping,  diagonal  damping  and  non- 
diagonal  damping  problems;  the  adjustment  of  the  basic 
formulas  to  accommodate  step  changes;  and,  step  size 
selection  criteria.  It  Is  shown  that  the  truncation 
error  concept  is  not  adequate  to  detect  numerical  In- 
stability for  general  structural  dynamics  analysis 
when  low  accuracy  is  requested.  The  "apparent  fre- 
quency" concept  is  Introduced  to  maximize  stable  step 
sizes  and  Is  shorn  to  be  stable  when  low  accuracy  Is 
requested. 

L.  INTRODUCTION 

The  central  difference  integration  formula  Is 
widely  used  for  direct  time  Integration  of  the  equa- 
tions of  motion  governing  structures.  This  popular- 
ity Is  due  to  its  favorable  fixed-step  stability  limit 
for  undamped  problems,  no  numerical  damping,  low  fre- 
quency distortion,  and  simple  implementation.  The 
formula  is  thus  well  suited  for  fixed-step  Integration 
cases  provided  the  step  size  chosen  does  not  exceed 
its  stability  limit  with  respect  to  the  highest  fre- 
quency of  the  system  under  consideration  and  the  so- 
lution accuracy  is  not  adversely  affected  by  such  a 
step  size  choice. 

In  practice,  however,  step  size  changes  are  often 
desired  to  select  a step  size  as  large  as  possible, 
consistent  with  a specified  local  error.  This  calls 
for  an  adaptation  of  the  central  difference  formula 
to  a variable-step  Integration  procedure.  There  are 
four  major  factors  that  affect  the  efficiency  of  a 
variable-step  Integration  procedure:  the  basic  Inte- 
gration formula  which  dictates  stability  and  accuracy 
for  fixed-step  Integration;  the  adjustment  of  the 
basic  formula  to  handle  atep  size  changes  and  to  min- 
imize any  loss  of  stability  and  accuracy;  the  step 
size  selection  strategy  to  satisfy  a specified  error 
bound  and  at  the  same  time  to  avoid  Instability;  and 
and  last,  but  not  least,  the  computer  Implementation. 

This  paper  Is  divided  Into  two  parts  to  address 
the  above  four  factors.  In  Part  I we  concentrate  on 
che  first  three  factors  as  they  are  closely  tied  to 
the  algorithmic  nature  of  the  Integration  process. 
Thus,  Part  I covers  theoretical  aspects  of  the  vari- 
able-step Integration  package.  Part  II  O)  describes 
the  computer  Implementation  and  the  performance  eval- 
uation of  the  package  through  carefully  chosen  numer- 
ical experiments.  The  many  "ifs"  to  treat  the  many 
"exceptions"  that  accompany  a typical  software  devel- 
opment process  are  also  presented. 


2.  OUTLINE  OF  PART  I 

The  stability  and  accuracy  characteristics  of  the 
fixed-step  central  difference  formula  are  well  under- 
stood for  undamped  equations  of  motion  governing  struc- 
tures, cf.,  Krelg  and  Key  (2)  or  Park  (3).  The  effect 
of  system  damping  on  Its  performance,  however.  Is  not 
well  understood  even  though  It  is  generally  acknow- 
ledged that  the  presence  of  nondiagonal  damping  re- 
duces Its  stability  margin.  In  Section  3 the  basic 
formulas  for  the  fixed-step  Integration  are  presented 
for  undamped,  diagonally  damped  and  nondlagonally 
damped  cases.  In  Section  4 the  analysis  of  stability 
and  accuracy  of  the  formulas  for  the  nondlagonally 
damped  equations  of  motion  Is  presented.  The  origin 
of  the  stability  reduction  due  to  nondlagonal  damping 
la  Identified  and  a means  to  minimize  this  stability 
reduction  Is  given. 

The  step  adjustment  techniques  are  described  In 
Section  5.  Two  techniques  are  introduced  to  adjust 
the  basic  Integration  formula  to  accommodate  step  size 
changes:  interpolation  of  pa3t  solutions  to  represent 
solution  vectors  at  equal  Intervals  of  the  present 
step  size  and  adoption  of  variable  coefficients  In  the 
integration  formula.  Care  Is  taxen  so  tnat  either 
technique  retains  the  characteristics,  as  close  as 
possible,  ot  the  fixed-step  formula. 

The  heart  of  the  present  paper  Is  concerned  with 
step  size  selection  strategies.  Several  variable-step 
Integration  procedures  have  been  proposed  for  explicit 
methods,  notably,  oy  Ktogh  (4),  Brayton  et  al.  (5), 
Gordon  and  Shampine  (b) , Zadunalsky  (V) , among  others. 
While  the  details  of  these  procedures  vary  consider- 
ably, their  common  presumption  Is  that  the  local 
truncation  error  Is  a reliable  measure  to  control  step 
size  changes.  Although  this  appears  to  be  the  case 
for  high-accuracy  requirements,  I.e.,  the  relative 
local  error  Is  less  than  10~5,  Its  applicability  to 
oscillatory  problems  such  as  structural  dynamics  prob- 
lems has  not  been  firmly  established.  Moreover,  In 
typical  structural  dynamics  analysis  much  lower  accu- 
racy (say,  a few  percent  relative  error)  Is  often  sat- 
isfactory due  to  analysis  cost  and  uncertainties  In 
structural  modeling.  In  Section  6,  the  applicable 
ranges  of  truncation  error-based  step  selection  cri- 
teria are  examined.  This  examination  Indicates  that 
for  low-accuracy  analysis  the  truncation  error  con- 
cept Is  unreliable.  The  step  size  Is  either  too  con- 
servative so  that  It  yields  more  accuracy  than  re- 
quested or  too  unconservative;  hence,  the  step  size 
Is  seldom  the  one  desired. 

In  Section  7 we  present  a new  atep  size  selection 
criterion.  The  criterion  Is  based  on  the  sampling 
concept  of  oscillating  signals,  In  which  frequencies 
are  determined  from  the  "apparent"  oscillations  of  the 
solution  components.  Here  we  have  been  guided  by 
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heuristic  schemes  ss  the  theory  Is  still  Incomplete. 
Nevertheless,  the  key  feature  of  the  new  criterion  Is 
shown  to  hsve  s sound  foundstlon  from  the  viewpoint  of 
the  Coursnt-Flsher  maximum-minimum  theorem  of  eigen- 
values (£) . the  numerical  experiments  presented  In 
Part  II  Illustrate  that  the  new  step  size  selection 
criterion  exhibits  fidelity  to  a prescribed  error 
bound  and  safeguards  against  numerical  Instability. 

A summary  of  the  new  results  presented  In  Part  I la 
given  In  Section  8. 

3.  BASIC  INTEGRATION  FORMULAS 

For  direct  time  Integration  of  the  discrete  equa- 
tions of  motion 


.n  1 ,.n+*j  , .n-4.  ... 

u ■ j (u  ♦ u ) (6) 

to  obtain  from  (4) 


.n+i 


n+1 


a1  Mn_i+h  [p-n  - 

n . . .n+i 


(7) 


where  £|  ■ £ + *1  h £ and  £2  " tl  * *1  h £.  Notice  that 
the  Inversion  of  the  matrix  Ej  presents  no  computa- 
tional difficult)  as  It  Involves  only  a combination 
of  two  diagonal  matrices.  The  stability  limit  of  (7) 
Is  known  to  be  the  same  as  the  undamped  case  (2)  • 


M u"  + £ u"  + 1 (u°)  - P (t") 


(1) 


governing  structural  dynamics,  the  summed  form  of  the 
central  difference  formula  for  fixed-step  integration 
Is 


.n+*i  .n-*j  , . -n 

u - u + h u 

n+1  n , . .n+Jj 

u - u + h u 


(2) 


where  u Is  the  displacement  vector,  h the  fixed-step 
size,  the  superscript  (_)n  indicates  the  vector  value 
at  the  discrete  time  tn,  the  superscript  dot  (J_)  tem- 
poral differentiation,  £ and  £ are  the  diagonal  mass 
and  damping  matrices,  l_  (u)  is  the  set  of  Internal 
forces  due  to  stiffness  that  oppose  the  structural 
motion  (stiffness  forces),  and  P (tn)  Is  the  applied 
load  at  time  tn.  For  a linear  problem  the  stiffness 
force  £ (u)  becomes 

I (u)  - K u (3) 


where  l(  is  the  linear  stiffness  matrix. 

A general  recursion  formula  Is  obtained  by  com- 
bining equations  (1)  and  (2) 


^n+i  . in* * + h jf 1 (P"  - £ u"  - f (u“) ) 


n+1  n , . .n+i 

u * u + h u 


(4) 


Equation  (4)  will  now  be  specialized  to  three  basic 
formulas  to  treat  undamped,  diagonally  damped,  and 
nondlagonally  damped  problems. 


Nondlagonally  Damped  Case 

In  this  case  formula  (7)  Is  not  applicable  as  It 
loses  the  simplicity  of  vector  operations  due  to  the 
attendant  nondiagonal  factoring  of  £j . In  order  to 
preserve  the  vector  calculation  advantage  of  the  ex- 
plicit formula  (2),  the  term  £ u"  needs  to  be  approx- 
imated. The  effects  of  approximating  £u  on  numeri- 
cal stability  has  been  evaluated  in  (2) . The  results 
show  that  as  the  damping  Increases  the  stability  limit 
of  the  central  difference  formula  Is  reduced  for  all 
cases  examined. 

The  possibility  of  minimizing  the  reduction  of 
the  nondlagonal  damping-induced  stability  margin  was 
Initially  studied  in  (9).  The  underlying  observation 
to  be  exploited  for  this  minimization  is  that  the 
stiffness  force  and  the  damping  terms  can  be  separa- 
ted, and  the  accuracy  of  the  stiffness  force  term  Is 
not  affected  by  the  choices  for  approximating  un. 

This  suggests  Iterations  on  the  damping  term  as  a 
profitable  avenue  if  the  gain  from  minimizing  the 
stability  reduction  outweighs  the  overhead  for  iterat- 
ing the  damping  term.  This  lain  general  the  case  for 
nonlinear  problems  as  the  evaluation  of  the  stiffness 
force  term  constitutes  the  bulk  of  the  computational 
effort  and  the  damping  term  often  remains  linear. 

In  principle,  iterations  on  the  damping  term  £ u 
will  improve  both  accuracy  and  stability^  In  practice, 
however,  the  number  of  Iterations  on  £ u has  to  be 
limited  because  the  increase  In  the  stability  margin 
and  accuracy  from  additional  iterations  has  been  shown 
to  diminish  beyond  the  first  two  Iterations  (9)  and 
there  Is  the  possibility  of  rejecting  the  solution  at 
that  very  time  step  even  If  the  damping  computation 
has- converged.  Therefore,  we  limit  the  Iterations  on 
the  damping  term  to  two  at  most  and  the  resulting 
basic  formula  takes  the  form  of 


Undamped  Case 

This  case  corresponds  to  £ ■ £ In  equation  (4) 
and  the  recursion  formula  becomes 


^n+i  . 4»-i+h„-l  [Pn-f  <„“)] 


n+1 


- u^hu1*1 


(5) 


It  Is  well  known  that  the 
joys  the  full  accuracy  of 
and  there  Is  no  reduction 
l"i(x  h i 2,  where  Is 
undamped  system. 


sbove  recursion  formula  en- 
the  integration  formula  (2) 
of  the  stability  margin  from 
the  highest  frequency  of  the 


Diagonally  Damped  Cass 

When  damping  Is  present,  difficulty  srlsss  In  the 
computet Ion  of  £ un  because  the  velocity  vectors  are 
computed  st  half -stop  points.  The  best  known  answer 
to  this  difficulty  Is  to  linearly  Interpolate  un  as 
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where  y and  g are  Integration  parameters  and  a is  an 
averaging  parameter  selected  to  minimise  the  loss  of 
the  stability  margin,  see  Stetter  (10) . 

We  now  examine  the  stability  and  the  accuracy  of 
formulas  (8)  through  (11)  for  the  nondiagonal ly  damped 
problem. 

4.  STABILITY  AND  ACCURACY  ANALYSIS  OF  FORMULAS  FOR 
NON-DIACONALLY  DAMPED  PROBLEMS 


The  analysis  procedure  for  stability  and  accuracy 
of  Integration  formulas  for  linear  (or  linearized) 
equations  of  motion  Is  well  known,  e.g.,  Gear  (11) 
and  Hall  and  Watt  (12) . For  structural  dynamic  equa- 
tions the  reader  Is  referred  to  Nickell  (13) . Krleg 
and  Key  (2).  Bathe  and  Wilson  (14) , and  Park  (15). 
among  others.  For  the  linear  case.  It  Is  well  known 
that  the  Integration  procedure  can  be  characterized 
by  considering  the  homogeneous  modal  form  of  equation 
(1)  as 


+ 2(ti>U  + Ul  u 


(12) 


where  ( is  the  damping  ratio. 

The  task  of  determining  suitable  sets  of  para- 
meters a,  6 and  y In  equations  (8)  to  (10)  was  pre- 
sented to  the  computer  (Bee  Appendix  A for  details); 
where  the  values  a - (0.0  ~ 1.0),  6 ■ (0.0  — 1.0)  and 
y - (0.0  ~ 0.5)  were  considered.  Of  all  combinations 
of  these  parameter  ranges,  four  sets  of  parameters, 
summarized  In  Table  1,  were  chosen  for  actual  imple- 
mentation due  to  their  favorable  stability  margins. 


Table  1 Candidate  formulas  for  nondlagonally  damped 
problems 


Formula 

Parameters 

FI 

a - 0.5,  6 - 1.0,  Y * 0 

F2 

a - 0.5  8 “ 1.0,  y - 0.5 

F3 

a * 0.75,  8 * 0.75,  y - 0.5 

F4 

a - 1.0,  8 * 1.0,  Y “ 0.0 

F5 

a • 0.0,  8 - 1.0,  Y • 0.0 

Figure  1 shows  the  stability  limits  of  the  four 
formulas,  along  with  the  conventional  formula  (F5), 

In  terms  of  the  normalized  step  size  uih  vs.  the  damp- 
ing ratio  £.  Observe  from  Figure  1 that  the  choice 
of  parameters  a,  8 and  y has  a dramatic  effect  on  the 
stability  limit.  For  example,  It  suggests  a smaller 
value  of  a In  equation  (10)  If  the  damping  ratio  (£) 
associated  with  the  highest  frequency  Is  greater  than 
0.5  and  vice  versa. 

Figure  2 shows  the  accuracy  characteristics  of 
formulas  (FI)  and  (F4)  in  Table  1 and  of  formula  (7) 
for  the  diagonally  damped  case  (see  Appendix  A for 
accuracy  definitions).  Notice  that  when  the  system 
la  lightly  damped,  l.e.,  C < 0.1,  formula  (F4)  is  more 
accurate  than  formula  (FI).  This  Is  because  within 
the  neighborhood  of  ( { 0.1,  formula  (F4)  enjoys  a 
larger  stability  margin  than  formula  (FI),  as  can  be 
observed  from  Figure  1.  Formulas  (F2)  and  (F3)  have 
more  or  less  similar  accuracy  characteristics  as  for- 
mula (FI)  as  their  stability  characteristics  are  quite 
similar.  The  accuracy  characteristics  of  Formula  (F5) 
are  not  Included  as  It  Is  the  least  desirable  from  a 
stability  viewpoint  (see  Figure  1).  It  Is  emphasized, 
however,  that  formula  (F5)  has  been  used  In  many  ex- 
isting explicit  Integration  packages  to  treat  non- 
diagonal  doping  In  the  equations  of  motion. 


Fig.  1 Stability  limits  of  basic  formulas 
(see  Table  1 for  formulas) 

From  the  foregoing  analysis  of  stability  and 
accuracy  of  the  four  basic  formulas  we  have  chosen 
formula  (FI)  for  heavily  damped  problems,  formula  (F4) 
for  lightly  damped  problons,  and  formula  (7)  for  diag- 
onally damped  problems.  In  practice,  however,  an 
accurate  estimate  of  the  damping  coefficient  associ- 
ated with  the  highest  frequency  of  the  system  Is  often 
difficult  or  even  Impossible.  Under  such  situations, 
formula  (FI)  Is  adopted  as  a default. 

So  far  In  this  section,  fixed-step  formulas  for 
various  system  characterizations  have  been  discussed. 
When  the  step  size  changes,  one  must  accomodate  the 
effect  of  the  step  size  change  on  each  basic  formula. 
This  will  be  dealt  with  in  the  next  section. 

5.  STEP  CHANGING  TECHNIQUES 

It  was  pointed  out  In  Section  2 that  step  size 
changes  can  be  accommodated  either  by  Interpolation 
of  past  solutions  or  the  use  of  vsrlable  coefficients 
In  the  Integration  formula.  In  general.  Interpolation 
requires  additional  solution  vectors  to  be  stored  In 
order  to  obtain  solution  vectors  at  equal  step  Inter- 
vals of  the  present  stap  size.  The  Introduction  of 
vsrlable  coefficients  alleviates  this  additional  stor- 
age requirement,  but  is  known  to  cause  Instability  if 
the  step  size  Is  changed  too  often  (16) . In  this 
study,  four  stap  changing  techniques  have  been  devel- 
oped and  they  are  summarized  In  Table  2.  Note  that 
step  size  changes  affect  the  stability  only  at  one 
step,  viz.,  the  time  step  at  which  the  stap  size  Is 
changed.  This  Is  because  the  formula  for  computing 
the  displacement  vector  Is  the  smae  for  all  step 
changing  techniques.  Thus,  in  order  to  examine  the 
effect  of  step  size  changes  on  the  local  stability  of 
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Fig.  2a  Accuracy  of  basic  formulas 
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Fig.  2b  Accuracy  of  basic  formulas 
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Fig.  2d  Accuracy  of  basic  formulas 


this  typa  of  Integration  procedure,  It  suffices  to 
examine  how  the  amplification  factor  X will  behave 
when  the  step  site  Is  changed.  The  reader  la  referred 
to  Appendix  I for  datalla  of  the  derivation  of  the 
local  stability  polynomiala,  and  the  results  are  sum- 
marised In  Figures  3 and  4 for  the  values  of  a*g  - 0.5 
and  y • 0 In  equations  (#)  through  (11).  Notice  from 
Figure  1 that,  whan  Increasing  the  stap  slsa,  tha 


velocity  averaging  technique  (FAVE)  la  most  desirable 
for  undamped  problems,  whereas  the  variable  coeffi- 
cient method  (FVCO)  outperforms  the  rest  for  heavily 
damped  problems.  When  reducing  the  step  site,  both 
the  extrapolation  (FEXT)  and  the  Interpolation  (PINT) 
techniquaa  become  ldertlcal  and  superior  to  the  FVCO 
for  lightly  damped  problems.  However,  the  egg-shaped 
unatable  tones  In  Figure  A for  heavily  damped  problem 


Table  2 Formulas  for  step  changes 


c 


Technique 

Formula 

Variable  Coefficient 
(FVCO) 

• nfj  .n-i  . 1,.  ,,  . ..n 

u - u + +h  ,)  u 

2 n n-1  — 

n+1  n .n+i 

u ■ u + h u 

n 

Average  Velocity 
(FAVE) 

un+i  - (3-r)  un-i  - (1-r)  un-*]  +hll“ 

l n 

n+1  * n .n+i 

u ■ u + h u 

n 

Extrapolation 

(FEXT) 

un+i  - if(7+r2)  4"“*  - (1-r2)  un-^)  + h iin 
o n 

n+l  n .n+i 

u B u + h u 1 

n 

Interpolation 

(FINT) 

un+*  - u"-J  + |(hn+hn_1)  [ (3+r)iin  + (l-r)un_1] 

n+l  n , . .n+i 

u ■ u + h u 

n 

1}  r " hn/hn-l 

2)  All  the  above 
when  r ■ 1. 

ormulas  recover  the  basic  formula  (2) 

DAMPING  RATIO  (() 

Fig.  3 Local  stability  of  various  step  Increasing 
techniques  (see  Table  2 for  formulas) 
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Fig.  4 Stability  of  various  step  reducing  tech 
nlques  (see  Table  2 for  formulas) 


suggest  that  one  should  use  the  FVCO  due  to  Its  small- 
er unstable  zone.  This  will  be  further  elaborated  in 
the  Evaluation  section  of  Part  II. 

6.  CONVENTIONAL  STEP  SIZE  SELECTION  STRATEGY 

In  this  section,  two  conventional  step  size  se- 
lection criteria  will  be  analyzed:  the  truncation 
error  and  the  difference  In  the  predictor-corrector 
solutions.  It  will  be  shown  that,  when  Jow  accuracy 
Is  requested,  these  criteria  are  either  unreliable 
and/or  too  expensive  to  use. 

Analysis  of  Truncation  Error 

In  order  to  make  our  point  clear,  we  will  use  the 
linear  homogeneous  undamped  problem,  viz., 

M un  + K u"  - 0 (13) 

which  can  be  decomposed  into  modal  components 

n _ 


1 + A £ 
where 
U ■ T £ 

i‘n  - L 
IT  *1  - A 


(14) 


n+1 


- 12  ”®x  {|ek+l1} 


(19) 


l<k<N 


Upon  multiplying  both  the  numerator  and  the  denominator 
by  u£+1,  equation  (18)  can  be  written  tn  terms  of  modal 
components  from  (14)  as 


N 2 2 n+1  A n+1 

n+1  . J-l  J 1 3 3 

ck  N n+1, 2 


X.  (tk 


n+1 


— '■u  9,  ) 

J-l  J j 

n+1  n 
q - q . 


(20) 


where  Aq 
Assuming  that 

l|Aatell2  i (2i) 

equation  (20)  can  be  approximated  in  the  form 

Kn+1  Kn 
n+1  . _K_  bk 


n+1 

ak 


where 


N 2 2 
Z (t.  q.T 

J-l  3 J J 


(22) 


(«Jl 


Z q,> 

J-l  J 3 


The  local  truncation  error  of  the  central  difference 
method  (2)  for  integrating  (13)  can  be  expressed  as 

n+1  h3  -n  h3  -1  „ .n 

- ■ 12 u ■ - n« 


h y—  I.,  f n+ 1 n . 

12  ^ " H > 


(15) 


Remark:  Assumption  (21)  Is  equivalent  to  the  approx- 
imation 


n+1 

ck 


,-n+l  -n, 

(UK  - V 

fr+1 


n+1  *n+l 
u.  u. 
k k 


, n+1, 2 
(uk  > 


, n,2 
(uk3 


(23) 


A proper  local  error  measurement  for  the  explicit  for- 
mula Is  to  evaluate  component -wise  relative  errors  and 
take  the  maximum  among  them,  viz., 


n+1 


||t1/u1|,  |t2/u2|,  ...  (16) 


It  must  be  noted  that  error  norms  of  ip-type 

lli/jillp  - {|ti/ui|p  + ...  + UIt/uNlp}1/p  <m 

although  they  are  suitable  for  implicit  formulas,  are 
not  applicable  to  explicit  formulas  because  they  tend 
to  mask  out  locally  unstable  components  which  propa- 
gate violencly  In  the  explicit  integration  process. 

To  simplify  the  subsequent  analysis  let  us  Introduce 
the  k-th  relative  error  component  e”**  as 


n+1 

ak 


•o  that 


component 

.••n+1  ..n, 

<UK  * V n+1  . „ 

k 


(18) 


ii {4 (a  + 4a)  ' I5* <a  + 4a)  1 (24) 

2 

We  now  define  the  k-th  "apparent  frequency"  (u)^)^  as 


<J)  - ^ 

• k \ 


(25) 


for  the  k-th  solution  component.  Then  It  follows  from 
the  Courant-Flscher  Mlnlmax  Theorem  (9)  that 


“min  4 (*4>k  4 W k ‘ 1 N 


(26) 


Equation  (22)  can  ba  Interpreted  aa  the  change  of  the 
apparent  frequency  (uij),  from  time  tn  to  tn+‘  for  the 


k-th  eolutlon  component,  via., 


c*1-  (Aw2) 


A k 


, 2,n-l  , 2,n 

<uA>k  * WA  k 


(27) 


ao  that  (19)  can  ba  axpreeaed  aa 
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n+1 


12  K^kLa,' 


1, 


(32)  - (34)  that  the  difference  dn**  between  the  pre- 
(28)  dieted  and  converged  solution  Is 


For  a fixed  error  bound  c,  equation  (28)  allows 

(h)2  < 12  e/(Au>2)  (29) 

max  A max  y 


while  the  maximum  step  size  (h  ) for  stability  is 

s max 


(h2)  < 4/u2 

s max  max 


(30) 


2 2 

For  stability  (h  )max  S (h8)max;  hence,  the  step  site 
selected  based  equation  (29)  must  satisfy 


2 2 

(Aw  ) > 3 e u) 

A max  max 


(31) 


Note  that  the  apparent  frequency  for  the  solu- 
tion component  u^  as  defined  by  equation  (25)  can  be 
viewed  as  a "component -wise"  Rayleigh  quotient.  This 
quotient  often  typifies  the  dominant  response  fre- 
quency for  u^  (not  necessarily  the  fundamental  one). 
Furthermore,  the  "component-wise"  Rayleigh  quotients, 
in  practice,  do  not  change  drastically  from  one  step 
to  the  next.  Consequently,  (6wj[)max  will  remain 
small.  From  (31)  observe  that  as  w2  increases, 
one  must  choose  a smaller  error  bound'  tor  stable  in- 
tegration. Therefore,  the  step  size  criteria  based 
on  the  truncation  error  concept  for  general-purpose 
structural  dynamics  analysis  Is  reliable  only  when 
high  accuracy  [small  e In  (29)]  Is  requested. 


Error  Analysis  of  a Predictor-Corrector  Pair 

The  central  difference  method  and  the  trapezoidal 
rule  will  be  used  as  a predictor-corrector  pair.  A 
widely  used  starting  procedure  is 


Predictor: 


u* 

~P 


.o 

u 


+ £h  u° 


u1 

V T 


+ h u 


(32) 


, n+1  n+1  n+1 

1 - ^ -up 


h 

4 


2 


...n+1  ..n, 

(u  - u ) 


(36) 


Note  that  this  difference  is  three  times  the  trunca- 
tion error  (15).  Therefore,  except  for  its  Increased 
err>  r constant  1/4  as  compared  to  1/12  In  equation 
(15),  the  analysis  and  conclusion  In  the  previous  sub- 
section on  "Analysis  of  Truncation  Error"  is  applica- 
ble to  equation  (36). 

The  foregoing  analysis  of  the  truncation  error- 
based  step  selection  strategy  and  the  predictor-cor- 
rector Iteration  procedure  present  structural  dynamics 
analysts  with  a real  dilemma  in  that  the  former  Is 
unreliable  and  the  latter  is  expensive  to  use.  This 
has  motivated  the  authors  to  develop  a new  step  selec- 
tion strategy,  described  in  the  next  section,  to  re- 
solve this  dilemma. 


7.  NEW  STEP  SIZE  SELECTION  STRATEGY 


It  has  been  demonstrated  In  the  previous  section 
that  the  use  of  the  conventional  step  size  selection 
criteria  based  on  the  truncation  error  concept  Is  lim- 
ited to  nonstiff  problems  or  to  a high-accuracy  Inte- 
gration requirement  (cf.,  more  than  five-digit  accu- 
racy). An  effective  use  of  explicit  methods,  thus 
calls  for  the  determination  of  the  highest  frequency 
of  the  system  In  order  to  limit  the  maximum  step  size. 
For  nonlinear  problems,  however,  this  often  presents 
a challenge  as  the  highest  frequency  can  change  signif- 
icantly. In  addition,  not  all  natural  frequencies  are 
excited  under  a given  load.  Therefore,  one  wishes  to 
compute  the  highest  apparent  (excited)  frequency  In 
order  to  maximize  the  stable  step  size.  To  this  end, 
a "perturbed  apparent  frequency"  w Is  Introduced 
In  the  form 


<"PA>1 


, n+1  ."n+1 
Au^  'Auj 

Au2 


(37) 


" if1  (1  - K M1)  (33) 

• o . h ...o  . ..1. 

■ u + j (u  + u ) 

(34) 

- u°  + h u°  + | (u°  + u1) 

The  solution  process  takes  the  pattern  of  prediction 
(P)  once  and  evaluation  (E)  and  correction  (C)  until 
the  solution  converges,  i.e.,  PlEC)*1  for  k-time  iter- 
ations. For  subsequent  steps,  (32a)  Is  replaced  by 
(2a) • k 

To  test  convergence  of  the  P(EC)  process  one 
usually  employs 


Evaluate: 


..1 


Corrector : 


for  the  i-th  component  of  the  solution  vector  u. 

For  the  linear,  homogeneous  undamped  problem  (13), 
the  "erturbed  apparent  frequency  (37)  can  be  expressed 
as 


. 2 , A 

‘“W  “ 


r 2 2 
. t 


r 2 .2 

A<,J 


li-Ji 


(38) 


where  (14)  Is  used  as  In  the  derivation  of  (20). 
Comparing  (38)  to  (25),  observe  that  the  perturbed 
apparent  frequency  (“pa)i  can  be  viewed  as  the  fre- 
quency associated  with  the  1-th  component-wise  noise. 
The  "maximum  perturbed  apparent  frequency"  is  now  de- 
fined as 


(35) 


where  C is  a convergence  threshold. 

Here  the  convergence  test  (35)  Is  used  for  ad- 
justing step  sizes.  Notice  that  the  central  differ- 
ence method  as  a predictor  provides  the  required  full 
accuracy,  and  the  trapezoidal  rule  as  a corrector  Is 
used  merely  to  check  stability.  It  can  be  shown  from 


(Vsa,  ' Mx  {Sa>1 SaV*  (39) 


Therefore,  (upA)max  Is  recognized  as  the  highest  com- 
ponent-wise noise  frequency  that  la  noticeable  (appar- 
ent); hence  («„,)  can  be  used  to  determine  h , 

» A max  max 


‘“Wmax 


(40) 
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for  stable  Integration. 

To  complete  the  step  size  selection,  (“p/^mln 
should  also  be  determined  so  that  the  number  of  sam- 
ples per  the  longest  period,  Ng,  satisfies  the  accu- 
racy requirement.  This  Is  easily  computed  from 

Ns  * ^‘Vsln"1  <*» 

Remark:  The  qualifier  "apparent"  Is  used  here  to  em- 
phasize only  those  frequencies  that  participate  In  the 
dynamic  response.  Note  (u>PA)max  la  the  highest  fre- 
quency among  the  active  frequencies  that  participate 
in  the  dynamic  response.  For  example,  a structure 
under  free-fall  gravity  motion  does  not  exhibit  any 
frequency  content  except  the  rigid  body  motion  (cf., 
zero  frequency).  In  this  case,  (wpA)max  zero  and 
the  step  size  determined  from  (wpA)max  muc^  larger 
than  the  critical  step  size  based  on  the  highest  fre- 
quency of  the  system  as  determined  from  the  eigenprob- 
lem. 

8.  CONCLUSIONS 

In  this  paper  the  term  "method"  is  used  to  In- 
clude the  integration  formulas,  the  manner  In  which 
the  step  size  changes  are  incorporated  Into  the  Inte- 
gration formulas  and  the  step  size  selection  strategy. 
We  now  summarize  the  major  results  of  the  present 
paper. 

1.  It  Is  expedient  to  evaluate  separately  the 
stiffness  force  term  and  the  damping  term.  For  heav- 
ily damped  problems,  two  evaluations  of  the  damping 
term  per  step  yield  a considerable  Improvement  on  the 
stability  margin  of  the  central  difference  formula 
(see  Table  1 for  details). 

2.  The  variable  coefficient  formulas  are  super- 
ior to  the  Interpolated  formulas  for  use  in  step 
change  adjustments  when  the  system  Is  heavily  damped. 
On  the  other  hand,  the  interpolated  formulas  become 
superior  when  the  system  is  lightly  damped  (see  Table 
2 for  details). 

3.  Step  size  selection  strategy  based  on  the 
truncation  error  concept  cannot  be  used  for  detecting 
the  critical  time  step  size  necessary  to  avoid  numer- 
ical Instability  when  low  accuracy  Is  requested. 

4.  The  proposed  "perturbed  apparent  frequency" 
criteria  to  govern  the  choice  of  the  maximum  possible 
step  size  at  each  step  has  been  shown  to  be  stable 
when  low  accuracy  is  requested. 
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APPENDIX  A 

CHARACTERISTIC  POLYNOMIAL  OF  INTEGRATION 
FORMULAS  FOR  THE  NONDIACONALLY  DAMPED  CASE 

Equations  (8)  through  (11)  can  be  expressed  in 
matrix  form  as 

, n n— 1 _ 

A J!  + ® X - o 

where 

n i-n  •n+*  n+li 
jr  - lu  , u , u J 

in  which 
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i::3 

[2£ug,  2 tug,  gu2l 
0 -10 
0 0 -ij 


(A.l) 


i - <*cn 


g - (1  - a)  y - § (1  - 8)  - o*Ytn 


0 • uh 

The  appropriate  elgenproblem  for  (A. la)  la  obtained 
by  taking 


n+1  n 

Z ■ 


(A. 2) 


which  la  aubatltuted  leto  (A.l)  to  yield 

(*£  - 1)  X.  * 0 (A. 3) 

Equation  (A. 3)  has  a nontrivial  aolutlon  only  If 

det  |14-1|  * 0 (A. 4) 

Expanding  (A. 4),  one  obtalna  the  following  character- 
lstlc  polynomial  equation 


A(A-l)2  + g R2!2  + 2tng  1(1-1) 
+ 2ttlg  (l-l)2  - 0 


(A.5) 


Stability:  In  order  for  (8)  to  give  a bounded  aolu- 
tlon, one  must  have 

lAl  3 1-0  (A. 6) 

which  Is  the  required  stability  condition. 

Accuracy:  Accuracy  can  be  assessed  by  comparing  the 
exact  solution 


. , U ± J X - c2)  Uh 


(A. 7) 


to  the  numerical  solution  associated  with  the  computed 
damping  ( 


J-  tc  ± j 


so  that  one  can  express 

frequency  error  (Z)  ■ 


(-3) 


x 100 


damping  error  (Z) 


(1  - «r/C>  x 100 


(A. 8) 


(A. 9) 


(A. 10) 


APPENDIX  a 

THE  EFFECT  OF  STEP  CHANCES  ON  LOCAL  STABILITY 
The  variable  coefficient  formula  from  Table  2 la 


■ u"  * ♦ J(h  ♦ h .)  u" 
n D”  1 


V1*1  - Un  ♦ h iT** 

n 


(8.1) 


which  after  eliminating  the  velocity  terms  results  la 

h2 

u"*1  - <1  ♦ r>  u"  ♦ r u"'1  - (1  + i)  «■“, 


Vhn-1 


vB.2) 


The  acceleration  u from  the  basic  formula  can  be  ex- 
pressed from  (11a)  for  the  modal  equation  (12)  In  the 
form 

-n  .2  k . _ . -n-1 

u - - g(u  u ♦ 2Cu  u ) - 2i>ug  u 

where 


g • 1 - aBou  h_/r 

n 


g ■ (l  - a)  1 - ^ - aB  you  h /r 

A H 


(B.3) 


For  Illustrative  purposes,  let  us  choose  v ■ 0,  6 ■ I 
ro  that  combined  difference  equation  of  (B.2)  and 
(8.3)  becomes 


u"*1  - (1+r)  (1-  g(j2  ♦ CHn)l  u" 


♦ (r  - (l-*T)g  {fi  ) u 
n 


n-1 


(8.4) 


Now,  we  use 


u"  - X u"-* 
c 


u"*1  - X u" 


XX  u 
c 


(8.5) 


where  X Is  the  amplification  factor  for  the  fixed- 
step  cal-<  as  described  In  Appendix  A.  However,  the 
amplification  factor,  X,  from  un  to  un+l  Is  different 
from  A . Substituting  (8.4)  Into  (8.4),  one  finds 


U+»>  11  - g(^  + tan), 


- (r  - (l+r)g  C 0 1/A 
n c 


(8.6) 


The  effect  of  step  changes  on  local  stability  can  now 
be  evaluated  from  (8.6)  by  requiring 

M s 

as  In  Appendix  A,  where  Ag  Is  determined  from  (A. 6). 
Other  cases  can  ba  similarly  evaluated  via  the 
procedure. 
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ABSTRACT 

The  varlable-atep  central  difference  method  de- 
veloped In  Part  I la  Implemented  aa  a atand-alone 
aoftware  package  that  la  eaally  acceaaed  by  exlatlng 
structural  dynamics  analyzers  (l.e.,  f lnl te-eleaent , 
-difference  discrete  element  computer  codea)  through 
a common  data  structure.  Input/output  (I/O)  manager 
and  a few  user-supplied  control  and  Interface  rou- 
tines. The  performance  of  this  package  la  evaluated 
through  a computer  study  of  four  sample  problems  that 
embody  a large  class  of  response  characteristics: 
linear  to  highly  nonlinear;  wave  to  structural  to  un- 
coupled component  response;  narrow  to  wide  frequency 
spread;  no  damping  to  over-damping  and  accuracy- 
critical  to  numerlcal-stabllity-crltlcal  response. 

The  present  method  la  successful  In  meeting  these 
response  characteristics  with  effectiveness  In  all 
the  above  cases  with  the  single  exception  of  pure 
wave  propagation. 

1 . INTRODUCTION 

The  development  of  a new  t’.ae  Integration  method 
for  structural  dynamics  analyses  not  only  requires  a 
theoretically  sound  method  but  also  requires  an  effi- 
cient implementation  and  evaluation  procedure  to 
determine  if  the  method  la  truly  viable.  In  this 
paper  the  Implementation  and  evaluation  procedure 
are  presented  with  special  emphasis  on  the  Interac- 
tion of  the  theoretical  development  with  the  Imple- 
mentation and  evaluation  procedure. 

The  variable-step  central  difference  method  des- 
cribed In  Part  I 0)  Is  Implemented  as  a stand-alone 
package  as  this  provides  the  greatest  possible  bene- 
fit to  users  of  the  method.  The  modular  design  of 
the  package  has  a minimal  effect  on  the  Interior  de- 
sign of  the  host  structural  analyzer  and  changes  In 
the  time-integration  method  are  Immediately  avail- 
able to  all  users. 

The  performance  of  a computational  method  Is 
very  dependent  on  the  Implementation  and  vice  verse. 
In  th'is  study  the  desired  performance  encompasses 
cost  efficiency,  reliability,  minimum  storage  and 
I/O,  and  easy  user  Interface.  These  performance 
goels  are  In  some  cases  at  odds  with  each  other, 
hence  several  compromises , necessery  to  achieve  the 
desired  perfoiaence,  were  made.  Fortunately,  these 
compromises  have  a minimal  effect  on  the  cost  effi- 
ciency and  reliability. 

To  evaluate  the  performance  of  the  Implementa- 
tion and  the  method  (the  two  cannot  be  separated) 
four  semple  problems  havs  been  generated.  The  sam- 
ple problems  have  a small  number  of  degreee  of  free- 
dom and  may  be  viewed  as  simple,  but  they  embody  a 


large  class  of  response  characteristic*  that  a suc- 
cessful method  should  be  able  to  treat  effectively. 
The  sample  problems  range  from  linear  to  highly  non- 
linear, from  wave  to  structural  to  uncoupled  compo- 
nent response,  from  narrow  to  wide  frequency  spread, 
from  no  damping  to  overdamping  and  from  accuracy- 
critical  to  numerlcal-stabllity-crltlcal.  With  the 
exception  of  the  pure  wave  propagation  all  the  char- 
acteristics are  effectively  treated  by  this  Integra- 
tion method.  For  one  of  the  sample  problems,  a model 
of  a drop  test  that  has  periods  of  free  fall  followed 
by  a rebound,  the  variable-step  method  achieved  a 
higher  average  step  size  than  could  be  used  for  a 
stable  fixed  step  calculation. 

The  paper  Is  divided  Into  Implementation  and 
evaluation  sections.  In  the  Implementation  section 
the  functional  form  of  the  stand-alone  package  Is 
discussed,  followed  by  a presentation  of  the  key 
details  Involved  In  Implementing  the  method  of  Part  I 
(1).  The  section  concludes  with  a brief  discussion 
on  performance  versus  Implementation.  The  numerical 
evaluation  section  presents  the  philosophy  of  the 
evaluation  procedure  followed  by  the  results  of  the 
evaluation.  In  this  section  many  of  the  decisions 
made  during  the  development  will  become  apparent  as 
the  effects  of  these  decisions  are  easily  seen. 

2.  IMPLEMENTATION 

In  this  section  the  basic  concept  and  Implemen- 
tation of  a stand-alone  time  Integration  package  and 
the  specific  Implementation  details  of  the  variable- 
step  central  difference  method  presented  In  Part  I 
are  presented.  The  emphasis  will  be  on  the  non- 
diagonal  damping,  the  error  evaluation  and  the  step- 
size  chauglng  strategy  techniques  as  these  are  the 
unique  features  of  this  method.  A brief  discussion 
of  performance  versus  Implementation  concludes  this 
section. 

Stand-Alone  Time  Integration  Package 

Figure  I presents  a schematic  overview  of  the 
modularized  software  that  has  evolved  during  the 
developawnt  of  structural  dynamic  time  integration 
procedures.  This  modularization  has  proven  to  be 
vary  cost  effective  In  that  existing  structural  dynam- 
ic analysers  are  easily  coupled  Into  the  time  Inte- 
grator software.  For  example,  various  techniques  to 
determine  the  nonlinear  terms  are  readily  evaluated 
by  Interfacing  with  structural  analyser*  using  dif- 
ferent nonlinear  evaluation  methods.  And  changes 
In  the  time  Integrator  are  Immediately  available  to 
all  .:he  structural  dynamic  analysers,  since  the  time 
ln.sgrator  Is  a stand-alone  package. 
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Fig.  1 Overview  of  structural  dynamic  response 
analysis 


The  term  "stand-alone  package"  means  the  Inte- 
grator Is  completely  self-contained;  It  Includes  the 
Integration  method  plus  a local  data  manager  that 
are  accessed  through  user-supplied  subroutines.  Indi- 
cated by  parentheses  In  Figure  1.  The  dashed  boxes 
in  Figure  1 Indicate  that  additional  interfaces  may 
be  required,  usually  to  provide  a change  in  the  data 
format  when  different  local  data  management  la  used 
by  the  structural  analyzer  and  the  time  Integrator. 

To  Interface  with  the  stand-alone  time  Integrator 
the  user  must  provide  five  subroutines  (some  may  be 
dummy  routines).  The  first  (DRIVER)  transmits  the 
problem  parameters,  the  diagonal  mass  matrix  {J.  (In 
the  case  of  diagonal  damping,  £ may  be  transmitted), 
and  the  initial  conditions  u°,  u°  to  the  time 
integrator.  Also  at  the  end  of  the  time  integration 
computations,  control  is  returned  to  DRIVER.  The 
second,  an  applied  loads  routine  (FORCE),  provides 
the  load  vector  data  to  the  time  integrator.  The 
third,  an  output  routine  (OUTPUT),  provides  display 
of  the  reaponse.  The  fourth  routine  f DFORCE)  pro- 
vides the  damping  force  as  a vector;  and  the  fifth 
routine  (SFORCE)  provides  the  stiffness  force  as  a 
vector.  Note;  All  data  is  transmitted  as  vectors, 
because  diagonal  matrices  are  treated  as  vectors. 

The  time  Integrator  treats  the  governing  struc- 
tural dynamics  equations  of  motion  at  the  n-th  dis- 
crete time  as 

Jjt  u"  d(6n)  + f (un)  - P(tn)  (1) 

with  the  Initial  conditions  u°  and  6°,  where 
M Is  the  diagonal  mass  matrix, 

7 is  the  sat  of  Internal  forces  due  to  damping  that 


oppose  the  structural  velocity, 
i_  la  the  set  of  Internal  forces  due  to  stiffness 
that  oppose  the  structural  displacement,  and 
P Is  the  applied  load. 

Note,  in  the  case  of  linear  diagonal  damping,  d - 
J)  un  (JO  a diagonal  matrix),  the  elements  of  £ (a 
vector)  may  be  transmitted  to  the  time  Integrator 
Instead  of  supplying  d.  The  vectors  d,  f_  and  P are 
the  ones  transmitted  to  the  time  integrator  through 
the  user-supplied  routines  DFORCE,  SFORCE,  and  FORCE. 
The  user  must  assure  the  efficient  computation  of 
these  vectors. 

Time  Integrator 

The  functional  form  of  the  stand-alone  time 
Integrator  package  developed  during  this  study  is 
Illustrated  In  Figure  2.  The  package  includes  ev- 
erything except  the  user-supplied  subroutines  which 
are  shown  to  provide  a reference  point  to  Figure  1. 
The  time  Integrator,  Inside  the  dashed  box.  Is  com- 
posed of  the  routine  CENDIF/VECTOR  and  the  associated 
subroutines  ROTATE,  ACLTN , ERROR  and  STEPSZ,  and  Is 
controlled  by  the  routine  STINT/CENDIF.  To  retain 
simplicity  the  time  Integrator  package  considered  In 
this  paper  Is  an  ln-core  operation,  but  It  can  easily 
be  extended  to  out-of-core  operation  for  very  large 
problems.  In  this  section  the  discussion  will  empha- 
size the  Items  In  the  time-loop  portion  of  the  time 
Integrator.  However,  some  other  Items  will  be  brief- 
ly considered  first. 

During  the  opening  phase,  the  initial  values 
specified  are  those  for  the  time  Integrator  logic, 
not  the  Initial  conditions  for  the  equations  of  mo- 
tion. These  values  Include  arrays,  counters,  and 
switches  for  fixed-  or  variable-step  Integration, 
for  printout,  etc.  The  starting  procedure  means 
taking  the  first  time  step  |l.e.,  correctly  starting 
the  integration;  see  Dahlqulst  (2)1.  The  Integrator 
starts  with  a step  ten  times  the  minimum  value  selec- 
ted by  the  user;  this  gives  the  Integrator  room  to 
cut  the  step  size  without  errorlng  off. 


i 


Fig.  2 Time  Integrator  functional  outline 
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A synopsis,  summary  Information,  is  Included  to 
provide  some  data  on  the  performance  of  the  time 
Integrator  for  the  problem  it  has  just  solved.  Print- 
ed at  the  end  of  each  run,  it  Includes  the  average 
time  step,  the  number  of  step  Increases,  the  number 
of  step  decreases,  the  number  of  time  steps,  the  num- 
ber of  evaluations  of  the  equilibrium  equations  and 
a distribution  count  of  the  step  sizes  used  during 
the  time  integration.  This  information  is  helpful 
in  determining  the  efficiency  of  the  calculatlona  so 
that,  as  experience  is  gained,  the  user  may  more 
effectively  choose  input  parameters  for  similar  prob- 
lems. 

The  time  loop  section  of  the  integrator  is  where 
all  the  work  is  done.  The  flow  of  computatlors  pro- 
ceeds from  Items  3 to  7 and  then  back  to  Item  3 until 
the  problem  is  solved  for  the  requested  time  range. 
Now,  each  item  in  the  time  loop  will  be  discussed. 

The  update  of  the  velocity  and  displacement  for 
each  degree  of  freedom  depends  on  the  conditions  pre- 
vailing at  that  degree  of  freedom  at  the  current 
step.  After  the  acceleration  at  time  tn  has  been 
successfully  computed  (see  below),  the  formulas  for 
computing  the  velocity  at  time  tn'*i  and  the  displace- 
ment at  time  tn+'  are  chosen  as  shown  in  Table  1. 

The  FVCO  and  FINT  formulas  were  chosen  for  their 
stability  characteristics  under  decreasing  and  in- 
creasing step  changes;  see  Part  I. 

From  Part  I it  is  seen  that  to  advance  the  step 
to  time  tn+1  the  acceleration  at  time  tn  must  be  de- 
termined. For  the  undamped  case  this  is  trivial; 
see  Equation  (5),  Part  I.  For  the  diagonal  damped 
case  a special  form  is  used  (see  equations  (6)  and 
(7),  Part  I]  and  the  computation  of  the  acceleration 


Table  1 Formulas  for  step  changes 


Condition 

Fornula 

1)  r • h /h  , i 1.0 

FVCO 

n n-1 

2)  uJ-0 

FVCO  for  conponant  J 

3)  |u"  - uJ'VluJl  > 0.01 

FVCO  for  component  J 

4)  r > 1.0 

FINT 

and  conditions  2 or  3 

do  not  apply 

la  not  directly  required.  But  to  evaluate  the  error 
the  acceleration  is  then  determined  by  using  equa- 
tions (7)  and  (8)  in  Part  I.  The  nondlagonally  damped 
case  is  treated  in  Part  I via  equations  (8)  - (11). 
This  implementation  is  based  on  > ■ 0 and  B “ 1 and 
a fixed  iteration  cycle  of  predict-correct-predlct. 

The  several  steps  shown  in  Part  I can  be  condensed. 

In  the  coding  two  arrays  dv  and  w are  used  in  con- 
junction with  a scratch  array,  workb,  in  the  follow- 
ing manner: 

M-‘  (Pn  - f (u")l 

<«  (rt 

6"'*  + $ ("n  ♦ Vj>  <dv  - ft  1 v>  <2> 


workb  • 

d(w) 

••n 

dv  - M * workb 
— «*»  ' 

u ■ 

where  a i»  as  defined  in  (10),  Part  I,  and  h and 
hn-i  -re  the  time  step  sizes  for  the  n-th  and 
(n-I)-th  steps.  In  (4)  the  Pn,  £(un)  and  d (x)  terns 
are  all  obtained  from  the  user  supplied  routines  and 
the  operation  in  the  third  of  (2)  with  w on  both 
sides  of  the  equal  sign  Is  the  standard  FORTRAN  re- 
placement operation. 

Now  that  the  acceleration  at  the  n-th  step  is 
known,  the  "perturbed  apparent  frequency"  measure 
can  be  applied  to  determine  if  the  solution  at  the 
n-th  step  is  acceptable.  The  error  is  computed  from 
equation  (37)  in  Part  I as 


• max(t", 


c ) 

MAXDEG  ’ 


where 


.2  n 
h a, 
n i 

4 »r 


(3) 


MAXDEC  is  the  size  of  the  solution  vector  in  equa- 
tion (1),  and  hn  is  the  n-th  time  step  size. 

If  the  denominator  in  (3)  Is  zero  or  if  |u?  - u9-*|/ 
lujl  < 1.6  x 10"8  for  |u"|  i 0.0,  the  error  is  set 
to  zero.  The  second  condition  eliminates  the  problem 
of  computing  an  error  based  on  machine  roundoff  in- 
stead of  system  response. 

Another  feature  that  Increases  the  computational 
efficiency,  but  requires  some  knowledge  of  the  prob- 
lem at  hand,  is  to  replace  the  denominator  term  b"  in 
(3)  by  1 


max 

1 1-J 


n- 1 


i) 


(4) 


where  nq,  is  the  (average)  bandwidth  of  the  1-th  deg- 
ree of  freedom.  In  the  present  study  mb  * 2 (MAXDEC/ 
10+1)  has  been  used  as  an  estimate  when  the  average 
bandwidth  was  not  easily  determined.  The  above  re- 
placement for  evaluating  the  error  (3)  may  reduce 
the  accuracy  of  some  high-frequency  components,  but 
it  Insures  that  they  do  not  grow  to  excess  and  it 
still  produces  an  accurate  efficient  solution  for 
the  dominant  response  modes. 

Once  the  error  measure  from  (3)  has  been  deter- 
mined a decision  must  be  made  as  to  what  step  size 
is  appropriate.  The  user  specifies  the  samples  per 
cycle  (N)  for  the  highest  apparent  frequency;  N must 
be  greeter  than  s,  as  s is  the  stability  limit  for 
the  central  difference  Integrator.  The  number  of 
samples  per  cycle  converts  to  the  error  measure  of 
(3)  as 


e • (»/N)2; 


(5) 


sae  Section  7 in  Part  I.  Based  on  the  ratio  of  the 
computed  error  (3)  to  the  desired  error  (5)  for  the 
highest  apparent  frequency,  l.e.. 
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P ■ £/«.  (6) 

the  step-size  selection  Is  made  as  shown  In  Table  2. 

In  Table  2 the  condition  for  Increase  must  be 
met  for  five  consecutive  steps  to  maintain  reliabil- 
ity. The  limits  on  the  step-size  change  are  neces- 
sary to  maintain  accuracy.  And  the  restriction  for 
S > 0.9  for  step  decreases  is  to  eliminate  the  possi- 
bility of  many  small  changes  at  one  time  step. 


Table  2 Step  size  selection  strategy 


Condition* 

Action 

0.1  < p < 1.0 

No  chania 

P < 0.1 

Incraaaa 

(for  tho  lut  5 
conaacutlva  stop*) 

new  old 

whara  S - win  l(*p)~^,  1.31 

p > 1.0 

Doers ass 

^naw  " **  N>ld 
whara  S - (5p)“1/5 

or  if  S < 2/3;  S • 2/3 

or  If  S > 0.9;  S - 0.9 

A step-size  decrease  Indicates  that  the  computed 
quantities  for  the  n-th  step  do  not  satisfy  stability 
or  the  accuracy  requirements  that  were  requested  by 
the  user, hence  this  step  is  rejected.  The  n-th  step 
Is  again  attempted  by  starting  at  Item  3 In  Figure  2; 
this  process  Is  continued  until  the  step  size  Is 
acceptable.  If  the  step  size  so  determined  Is  less 
than  the  minimum  value  set  by  the  user,  the  code  re- 
turns control  to  the  user.  If  the  step  size  Is  de- 
creased more  than  four  times  at  one  time  step,  a 
restart  of  the  problem  Is  attempted  at  Item  2 with 
the  last  accepted  displacement  and  velocity  vectors 
used  as  Initial  conditions.  If  this  falls,  control 
returns  to  the  user.  At  this  point,  the  user  must 
reconsider  the  control  parameters  and  begin  again. 

If  the  step  size  Is  acceptable,  the  displace- 
ment and  velocity  are  transmitted  to  the  user  for 
display  purposes  through  the  routine  OUTPUT.  Because 
the  velocity  Is  computed  at  half  steps,  extrapolated 
values  at  whole  steps  are  provided  via  the  formula: 

u"  • 6n'S  + S h u"  . (7) 

_ — n — 

where  h„  Is  the  n-th  step  size. 

It  should  be  aaiphaslzed  that  many  of  the  deci- 
sions Involved  In  this  Implementation  are  based  on 
maserlcal  experiments  with  a limited  class  of  prob- 
lesis.  Hence,  for  another  set  of  problems,  one  might 
utilise  a different  set  of  decisions.  The  subjective 
decisions  are:  1)  the  coupling  concept  to  determine 
the  denominator  In  equation  (3);  2)  the  limits  on  p, 
especially  the  lower  one;  and  3)  the  factors  1.3, 

2/3  and  0.9  that  appear  In  Table  2. 

The  coding  required  for  the  varlable-atep  central 
difference  time  Integrator  la  obviously  more  Involved 
than  that  for  a simple  fixed  step  version.  But  the 
cost  differences  In  computing  with  a fixed  step  ver- 


sus variable  steps  are  minimal,  especially  when  one 
considers  the  power  of  the  variable  step  method.  In 
some  esses  the  variable-step  method  has  been  found 
to  be  cheaper  than  the  fixed-step  method  even  If  the 
maximum  fixed-step  size  could  be  accurately  determined 
a priori. 

Performance  versus  Implementation 

The  performance  of  a time  Integration  software 
package  can  be  considered  from  two  perspectives:  1) 
the  overall  performance  of  the  complete  package;  and 
2)  the  performance  of  the  time  Integrator.  The  first 
addresses  the  ease  of  use,  I/O  requirements  and  over- 
all efficiency,  whereas  the  second  Is  related  to  the 
computational  efficiency  of  the  method.  If  the  per- 
formance Is  poor  in  either  area,  the  usefulness  of 
the  package  Is  weakened. 

In  the  authors'  programming  environment  the 
stand-alone  package  Implementation  has  greatly  In- 
creased performance.  Usually,  new  problems  require 
the  adaptation  of  existing  software  to  produce  some 
new  analysis  capability.  With  an  existing  opera- 
tional time  Integration  package  and  using  a coamon 
I/O  manager  for  all  software,  a transient  response 
analysis  package  can  be  built  very  quickly  by  very 
simply  plugging  modules  together.  Of  course.  In  a 
strictly  production  analysis  environment,  a time  In- 
tegrator embedded  In  the  coding  may  be  more  desirable; 
see  Fellppa  and  Park  (3)  for  further  discussion  of 
this  subject. 

For  the  time  integrator  the  performance  Is  a 
strong  function  of  the  reliability  of  the  method.  A 
totally  reliable  variable-step  time  Integrator  would 
be  costly,  but  In  a majority  of  cases  near  total  re- 
liability with  modest  costs  can  be  obtained.  For 
this  Implementation  a large  number  of  samples  per 
cycle,  approximately  4s,  plus  error  evaluation  based 
on  each  degree  of  freedom  will  produce  close  to  total 
reliability.  The  most  efficient  (time  and  cost  of 
computer  operation)  Is  obtained  with  a sampling  of 
* and  error  evaluation  based  on  coupled  degrees  of 
freedom.  Experience  has  shown  that  this  later  ap- 
proach Is  roughly  802  reliable,  l.e.,  the  integration 
remains  stable  and  produces  a reasonably  accurate 
solution  (101  error  at  most). 

In  the  coding  of  the  method,  the  biggest  compro- 
mise wss  made  In  core  storage  used  for  the  computa- 
tions. To  keep  the  cost  doMi,  It  wss  found  that 
using  scratch  arrays  to  cut  down  on  the  number  of  DO- 
loops,  and  saving  the  three  most  current  displace- 
ments, velocity  and  acceleration  vectors  to  enable 
backtracking  without  restarting  and  for  error  eval- 
uation was  requlrad.  This  limits  the  problem  size 
that  can  be  solved  In-core  to  between  700  and  1800 
degrees  of  freedom  depending  on  the  damping  option 
and  time  step  option.  This  Is  still  a fairly  respect- 
ably-slzed  problem  for  dynamic  analysis,  so  it  la  not 
considered  limiting. 

3.  EVALUATION 

The  development  of  a new  Integration  method  and 
its  computer  Implementation  nuat  be  complimented  by 
a series  of  evaluation  tests  to  assess  Its  effec- 
tiveness In  engineering  analysis  environments.  Four 
simple  structural  dynamic  models  have  been  chosen 
that  cover  a wide  range  of  dynamic  characteristics 
that  are  useful  In  evaluating  various  algorithms  and 
Implementations.  In  this  section  the  performance  of 
the  variable-step  central  difference  method  will  be 
Illustrated  using  these  models. 

The  four  models  are  a four-wheel  dolly,  a drop 
taat,  a cantilever  beam,  and  an  axial  bar;  these 
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are  described  In  detail  in  Appendix  A.  The  four- 
wheel  dolly  model  has  a relatively  narrow  band  of 
frequency  components,  roughly  201  damping  and  the 
response  Is  characterized  by  an  Initial  transient, 
followed  by  a damped  free  vibration.  The  drop  test 
model  also  has  a narrow  band  of  frequency  components, 
102  and  1102  damping  elements  and  the  response  Is  the 
highly  nonlinear  motion  of  a package  falling  to  the 
ground.  Impacting  the  ground,  rebounding,  falling. 
Impacting,  etc.  The  cantilever  beam  Is  a simple  beam 
bending  model  Including  translational  and  rotational 
degrees  of  freedom,  hence  a fairly  stiff  (wide  fre- 
quency range)  problem.  Both  linear  and  nonlinear 
bending  characteristics  are  considered  for  the  un- 
damped response  to  a tip  step  load.  The  axial  bar  Is 
a simple  model  of  wave  propagation  in  a free-free  bar 
subjected  to  an  initial  displacement.  These  models 
cover  a wide  range  of  characteristics  encountered  In  a 
structural  dynamic  analysis:  linear  to  highly  non- 
linear; wave  to  structural  to  uncoupled  component  re- 
sponse; narrow  to  wide  frequency  spread;  no  damping  to 
overdamping;  and  accuracy-critical  to  numerlcal-sta- 
blllty-crltical  response. 

i'-mi-  Wheel  Dolly  (Moderate  Damping,  Narrow  Frequency 
Spread) 

The  computed  displacements  of  wheels  1 and  4 and 
the  force  between  ground  and  wheels  1 snd  4 for  the 
four-wheel  dolly  are  shown  In  Figure  3.  Note  the  dis- 
placements exhibit  early  time  high  frequency  and  late 
time  low  frequency.  The  spring  force  at  wheel  1,  FK5, 
shows  that  wheel  1 lifted  off  the  ground  between  0.02 
and  0.04  seconds  and  again  for  an  Instant  Just  after 
0.1  seconds  and  the  spring  force  at  wheel  4,  FKg,  shows 
that  wheel  4 lifted  off  the  ground  between  0.07  and 
0.09  seconds. 

The  time  step  size  histories  for  different  accu- 
racy specifications  are  shown  In  Figure  4.  The  highest 
average  time  step  was  obtained  with  a sampling  of  x 
samples/cycle  of  the  highest  apparent  frequency;  this 
Is  the  minimum  sampling  to  obtain  numerical  stability. 
The  lowest  average  time  step  was  obtained  with  a sam- 
pling of  3*  samples/cycle  of  the  highest  apparent 
frequency;  this  provides  a very  accurate  solution. 

Both  computations  were  made  with  a • 1.0;  see  (2)  and 
the  discussion  below.  The  middle  average  time  step 
size  shown  In  Figure  4 was  obtained  from  a truncation 
error  control  method;  see  Part  1 (1).  Note  that  both 
3x  sampling  apparent  frequency  and  truncation  error 
methods  used  a relatively  small  step  size  and  actually 
reduced  the  step  size  some  during  the  Initial  transient 
nonlinear  response.  This  Indicates  that  accuracy  is 
controlling  the  Integration  during  this  phase.  The 
3*  sampling  computations  show  that  after  the  Initial 
transient  the  step  size  variation  settles  down  but 
does  not  Increase  very  much  as  the  requested  accuracy 
is  controlling  the  step  size.  The  it  sampling  computa- 
tions Illustrate  that  the  step  size  generally  increas- 
es throughout.  The  step  size  Is  Increased  very  rapid- 
ly and  maintains  a large  step  size  resulting  in  a very 
quick  computation  that  Is  still  quite  accurate  (within 
52  of  the  most  accurate  computation  made  for  this  prob- 
lem). These  curves.  Figure  4;  show  the  effectiveness 
of  the  perturbed  apparent  frequent  measure.  On  the 
other  hand,  for  a low-accuracy  but  stable  run  with  the 
truncation  error  concept,  the  best  obtainable  is  the 
one  shown  in  Figure  4.  If  a lower  accuracy  were  re- 
quested it  would  increase  the  step  else  quickly  but 
would  go  unstable  at  later  tlstts;  if  a higher  accuracy 
were  requested  the  average  step  site  would  be  painfully 
saull.  Also,  the  truncation  error  method  doea  not  pro- 
duce the  consistency  that  is  shown  by  the  apparent 
frequency  method. 
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Fig.  3 Four-wheel-dolly  response 
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Fig.  4 Four-wheel-dolly  step  size  histories 

In  Part  I the  Influence  of  the  parameter  a on 
stability  and  accuracy  was  analyzed;  this  Influence 
is  verified  In  Table  3,  which  summarizes  some  results 
for  the  four-wheel-dolly  response.  The  effect  of  var- 
ious values  of  a on  the  average  step  size,  computer 
time,  and  the  accuracy  of  the  maximum  displacement  of 
wheel  l and  the  minimum  acceleration  of  wheel  2 (cho- 
sen as  they  are  the  most  sensitive)  is  considered. 

For  both  sampling  rates  (it  and  3x)  the  average  step 
size  and  hence  the  computer  time  is  essentially  the 
same  indicating  that  a does  not  Influence  the  stabil- 
ity limit  for  a moderately  damped  (approximately  202) 
problem;  this  la  in  agreement  with  Figure  1 of  Part  I. 
On  the  other  hand,  the  accuracy  is  Influenced  by  the 
choice  of  a;  the  values  given  for  displacement  and 
acceleration  for  N ■ 3ir,  a ■ 1.0  are  the  moat  accurata 
and  compare  very  well  with  other  known  accurate  inte- 
gration methods.  Notice  that  in  all  cases  the  accu- 
racy Improves  as  a Increases;  this  la  in  agreement  with 
Figures  2a  and  2c  of  Part  I.  This  problem  was  also 
run  for  two  fixed  atep  cases  corresponding  to  the 
average  time  atep  for  the  N • » and  3«,  a - 1.0  cases. 
It  is  seen  that  the  fixed  step  computations  are  aesen- 
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Table  3 Four-wheel-dolly  summary 


( 


Saaple 

• 

At 

avaraga 

(aac) 

CAU 

(aac) 

TOTAL 

(aac) 

(cm) 

(•/..cJ) 

* 

0.4 

0.0049261 

1.111 

11.996 

1.81721 

-2.11092 

• 

0.1 

0.0049261 

1.112 

11.910 

1.81919 

-2.41147 

• 

0.6 

0.0046948 

1.489 

16.109 

1.81042 

-2. 10810 

a 

0.8 

0.0049261 

1.140 

11.912 

1.81704 

-2.11868 

a 

1.0 

0.0048780 

1.  161 

11.912 

1.8641o 

-1.89901 

)i 

0.4 

0.0020619 

6.280 

19.981 

1.81119 

-2.27219 

la 

0.1 

0.0020661 

6.280 

19.947 

1.81718 

-2.17612 

la 

0.6 

0.0022026 

1.910 

19.484 

1.81911 

-2.11606 

la 

0.8 

0.0021101 

6.084 

19.671 

1.86161 

-2.01129 

la 

1.0 

0.0020661 

6.270 

19.861 

1.86807 

-1.99142 

Flaad  Stap 

1.0 

0.0048780 

1.171 

16.118 

1.91211 

-1.91671 

Flaad  Stap 

... 

0.0020661 

6.788 

21.192 

1.87219 

-1.98771 

daily  equivalent  in  computer  time  to  the  variable 
step  computations;  hence  there  Is  minimal  overhead  In- 
volved In  using  the  variable  step  feature.  Secondly, 
note  that  the  accuracy  for  the  larger  step  size  run 
la  poorer  than  for  the  comparable  variable  step  run. 
This  occurs  because  the  variable  step  run  used  a smal- 
ler Initial  time  step  when  the  response  is  accuracy 
critical.  There  appears  to  be  little  reason  to  recom- 
mend a fixed  step  computation. 

Drop  Test  (Heavily  Damped,  Component  Responses) 

Two  representative  response  histories  for  the 
drop  test  problem  are  shown  In  Figure  5;  the  dis- 
placement of  degree  of  freedom  7,  the  base,  and  the 
acceleration  of  degree  of  freedom  5,  the  10Z  damped- 
100  Hz  oscillator.  When  the  displacement  is  less 
than  0.0  the  base  Is  in  contact  with  the  nonlinear 
spring  representing  the  ground,  the  positive  slope 
portions  of  the  displacement  response  Indicate  rebound 
and  the  negative  slope  portions  Indicate  free  fall  con- 
ditions for  displacements  greater  than  0.0.  The  res- 
ponse Is  like  a bouncing  ball.  Whenever  the  structure 
contacts  the  ground  a shock  la  felt  by  the  oscillators. 
Indicated  by  the  sharp  spikes  seen  on  the  acceleration 
response.  These  spikes  dampen  out  and  during  the  free 
fall  a constant  free  fall  acceleration  of  -9.8  m/s^  Is 
obtained. 

The  step  size  history  for  a calculation  with  3* 
samples/cycle  of  the  higheat  apparent  frequency  and 
a ■ 0.25  la  shown  In  Figure  6.  During  the  Initial 
free  fall  the  step  size  Increases  very  rapidly,  since 
the  apparent  frequency  of  thia  Initial  free  fall  la 
zero.  Once  the  Initial  Impact  occurs  the  step  size  la 
decraaaad  and  after  this  the  step  size  never  returns 
to  Its  early  high  value  since  the  oscillations  of  tha 
various  components  restrict  the  maxima  step.  It  la 
generally  seen  that  the  step  site  la  decreased  upon 
Impact  and  lncreasas  during  free  fall,  although  the 
general  pattern  la  somewhat  clouded  by  the  very  diffi- 
cult accuracy  and  stability  conditions  presented  by 
the  overdamped  degrees  of  fraedom. 

This  problem  can  be  successfully  run  for  a flxad 
time  stap  of  0.001  sec,  but  for  0.0013  sec  the  calcula- 
tions become  unstable.  From  Figure  6 It  la  seen  that 
the  minimus  atep  slza  la  slightly  laas  than  0.001  aac 
and  this  appaars  to  be  dictated  by  stability  consldar- 
atlona  baaed  on  the  fixed-step  calculations.  Hare  the 
varlab 1 e-atep  method  la  clearly  superior  to  the  fixed- 
step  method.  The  variable-step  method  takea  advantage 
of  the  Instantaneous  conditions  to  run  with  a near 


Fig.  5 Drop  test  response 
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Fig.  6 Drop  test  step  size  history 


optimum  step  size  while  the  fixed  step  size  must  be 
restrlctlvely  small  to  maintain  stability  during  crit- 
ical portions  of  the  response. 

In  the  preceding  section,  "Time  Integrator",  the 
use  of  the  FINT  formula  for  Increasing  step  sizes 
under  constant  or  near-constant  acceleration  was  pre- 
sented. The  need  for  this  computational  detail  was 
made  apparent  during  the  evaluation  studies  of  the 
drop  test  problem.  During  the  Initial  free  fall  under 
constant  gravity  acceleration,  the  calculations  exhib- 
ited a mild  Instability  which  contaminated  the  subse- 
quent calculations  and  error  control  when  the  FVCO 
formula  was  used  to  compute  the  response.  By  switch- 
ing to  the  FINT  formula  to  take  advantage  of  the  more 
uniform  stability  limit  for  all  damping  ratios  (see 
Figure  3 of  Part  I)  this  problem  area  was  eliminated. 

The  behavior  of  the  drop  test  model  response  la 
summarized  in  Table  4 for  various  samplings  and  a's. 

For  a fixed  a (a  ■ 0.25)  the  average  step  size  de- 
creases for  Increased  sampling  rate  In  a very  consis- 
tent manner.  For  the  2s  sampling  rate,  variations  In 
a cause  large  variations  In  the  maximum  computed  accel- 
eration. For  the  lower  values  of  a the  maximum  accel- 
eration occurs  at  degree  of  freedom  5,  the  high-fre- 
quency lightly-damped  oscillator.  But  for  the  higher 
values  of  a the  maximum  acceleration  occurs  at  degree 
of  freedom  6,  the  high-frequency  overdamped  oscillator. 
This  Indicates  that  for  heavily  damped  problems  a fair- 
ly small  a should  be  uaed  to  preserve  accuracy.  This 
result  Is  foreseen  In  Part  I,  but  from  the  maurlcal 
experiments  It  Is  seen  that  an  a samller  than  that  con- 
sidered In  Part  I should  be  used.  The  acceleration 
results  show  that  for  a low  sampling  rata  or  a larger 
a the  overdamped  reaponae  la  tha  moat  difficult  to 
control  from  both  stability  and  accuracy  vlawpolnta. 
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Table  4 Drop  test  summary 


Staple 

a 

At 

(hc) 

CAU 

(••c) 

TOTAL 

(••C) 

|X| 

'mi 

<C«) 

1*1— 

(■7m'| 

• 

0.24 

0.001844) 

4.741 

18.887 

4 ). 208* 

1092.911* 

J*/2 

0.24 

0.0017889 

4.912 

19.0)4 

43.172. 

)79.6795 

2« 

0.24 

0.0017442 

6.074 

19. 313 

4 3.071 j 

164 . 338^ 

2i 

0. 10 

0.0019194 

4.487 

18.778 

43.007J 

139.6825 

2* 

0.40 

0.0024096 

4.492 

17.3)7 

42.934J 

279. 

2» 

0.40 

0.002)041 

4.772 

17.443 

42. 890 j 

1023. 4936 

In 

0.24 

0.0017094 

6.189 

16.91) 

42. 949 j 

123. 400^ 

4n 

0.24 

0.0016779 

6.197 

19.478 

42.9)4 j 

1 16.034^ 

• Subscripts  Indicate  degree  of  f reedoa 


The  displacement  results  show  that  for  a low-frequency, 
lightly-damped  element  very  good  accuracy  Is  obtained 
for  any  combination  of  sampling  and  o.  The  major  con- 
clusion Is  that  for  heavily  damped  problems  an  i > 

0.25  ■*  0.30  should  be  used  and  a sampling  rate  of  at 
least  2n  should  be  used  to  maintain  accuracy  and  sta- 
bility. This  Is  especially  true  when  the  highest 
frequency  component  is  also  one  of  the  heavily  damped 
components  as  In  this  example. 

Cantilever  Beam  (Wide  Frequency  Spread) 

The  tip  rotation  and  displacement  (degrees  of 
freedom  1 and  2)  for  the  linear  cantilever  beam  re- 
sponse Is  Illustrated  In  Figure  7 for  variable-step 
computations  based  on  the  apparent  frequency  and  the 
truncation  error  step  control  methods.  For  all  sam- 
pling rates  of  » or  greater,  the  response  computations 
based  on  the  apparent  frequency  method  are  essentially 
Identical  for  displacements;  only  the  accelerations 
exhibit  a converging  pattern  for  increasing  sampling 
rates.  This  Is  due  to  the  widely  spread  frequency 
content  of  this  problem  in  which  the  overall  displace- 
ment response  Is  composed  of  low-frequency  components 
while  the  high-frequency  components  are  controlling 
the  step  size.  Figure  7 shows  the  truncation  error 
control  rotational  response  has  a very  high  frequency 
hash  of  relatively  large  magnitude  superimposed  on 
the  actual  response.  This  hash  appears  because  the 
truncation  error  control  was  unable  to  detect  the 
stability  limit  until  the  step  size  exceeded  the  sta- 
ble step  size.  It  was  found  that  for  any  specified 
error  of  reasonable  magnitude  the  truncation  error 
control  method  was  unable  to  maintain  a stable  Inte- 
gration for  the  cantilever  beam.  For  extremely  low 
error  (high  accuracy)  the  truncation  error  control 
remains  stable  but  the  average  step  size  Is  at  least 
an  order  of  magnitude  less  than  the  stable  step  size. 

On  the  other  hand,  for  the  s sampling  rate  the  appar- 
ent frequency  method  achieves  an  average  step  size 
that  is  roughly  85X  of  the  stability  limit.  Also 
this  method  Is  able  to  Increase  the  time  step  very 
rapidly  up  to  the  maximum  step  size  used  and  maintain 
this  step  size  for  the  majority  of  the  computations. 
Hence,  the  average  step  size  Is  nearly  equal  to  the 
maximum  step  size. 

The  results  for  various  sample  rates  for  the  lin- 
ear and  nonlinear  cantilever  beam  response  are  summar- 
ized In  Table  5.  Notice  that  as  the  sampling  rate  In- 
creases the  average  step  size  decreases  in  a consis- 
tent manner.  The  displacement  la  converged  for  all 
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Fig.  7 Cantilever  beam  tip  rotation  and  displace- 
ment response 


Table  5 Cantilever  beam  summary 


Staple 

At 

(W««c) 

CAU 

(etc) 

TOTAL 

(eec) 

M 

(ca) 

m 

(ka/eec2) 

Linear  C 

ntl lever 

eaa 

If 

0.66445 

7.802 

21.266 

0.08878 

6676.64 

3 »/Z 

0.47962 

10.176 

24.093 

0.08877 

7998.71 

2" 

0.33113 

13.775 

28.410 

0.08877 

7131.56 

3ir 

0.29445 

15.196 

30.162 

0.08877 

7224.52 

4* 

0.18657 

22.840 

39.794 

0.08877 

6694.93 

Nonlinear 

Cantilever 

Beam 

w 

0.61728 

9.349 

23.214 

0.11575 

7503.16 

3»/2 

0.51414 

10.664 

24.827 

0.11575 

5304.03 

2 9 

0.36496 

13.941 

28.768 

0.11573 

5688.08 

sampling  rates  and  the  acceleration  shows  some  ran- 
domness but  the  deviation  la  very  small  especially  if 
the  comparison  is  made  among  the  higher  sampling  rates. 
The  nonllnesr  response  shows  that  softening  of  the 
bending  properties  of  the  beam  generally  Increases  the 
average  atep  size  as  would  be  expected.  The  only  dis- 
crepancy In  theae  observatlona  la  seen  at  the  * aam- 
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CONCLUSIONS 


pllng;  this  Is  probably  due  to  slight  steps  over  the 
stability  boundary.  The  parameter  a has  no  Influence 
on  the  undamped  problem,  so  it  is  not  considered. 

Axial  Bar  (Wave  Propagation) 

The  exact  solution  for  an  initial  displacement  of 
2.54  cm  at  the  center  node  of  the  bar  (actually  a tri- 
angular spatial  displacement  of  0.0  at  the  nodes  on 
either  side  of  the  center  and  2.54  cm  at  the  center) 
is  two  triangular  displacement  pulses  of  half  the  orig- 
inal magnitude,  one  traveling  in  the  negative  direction 
and  the  other  in  the  positive  direction.  At  the  free 
boundary  this  pulse  is  reflected.  If  this  model  is 
run  with  the  time  step  equal  to  the  critical  time  step 
(fixed-step  computation)  the  numerical  results  are 
exact,  as  shown  in  the  left  half  of  Figure  8 for  the 
displacement  response  at  nodes  11,  14,  and  17.  Note 
the  magnitudes  and  arrival  times  are  exact.  The  varia- 
ble-step method,  presented  In  this  paper,  produces  the 
correct  results  if  one  starts  with  the  critical  time 
step  and  uses  a sample  rate  of  n;  l.e.,  the  time  step 
is  never  Increased  or  decreased.  Hence,  the  apparent 
frequency  concept  functions  properly  for  the  wave 
propagation  problem.  But  at  any  other  initial  time 
step  the  error  computation  detects  an  error  and  then 
the  step  size  is  always  decreased.  The  response  is 
the  hashy  Junk  shown  in  the  right-hand  half  of  Fig- 
ure 8.  The  problem  of  the  proper  error  control  for  a 
pure  wave  propagation  response  has  not  been  resolved. 
This  will  have  to  consider  the  shape  of  the  error- 
step  size  curve  so  that  the  slope  of  the  curve  can 
enter  into  the  decision  whether  to  Increase  or  de- 
crease the  step  size  to  reduce  the  error. 
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Fig.  8 Axial  bar  response 


The  implementation  of  the  variable-step  central 
difference  method  developed  in  Part  I as  a stand-alone 
software  package  has  been  shown  to  be  a viable  product. 
The  package  is  cost  efficient  and  easily  accessed  by 
existing  structural  analyzers.  The  evaluation  of  the 
method  showed  that  the  apparent  frequency  concept  is 
vastly  superior  to  the  truncation  error  for  control- 
ling stability  and  accuracy.  Of  the  wide  range  of 
problem  characteristics  considered  only  pure  wave  prop- 
agation Is  not  properly  controlled  by  the  step  size 
selection  strategy  presented  here.  The  wave  propaga- 
tion control  problem  is  being  investigated  and  the 
resolution  of  the  problem  will  be  reported  in  a future 
publication. 
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APPENDIX  A 
EVALUATION  MODELS 

The  discrete  models  used  to  evaluate  the  time 
Integrator  are  presented  in  this  appendix.  A brlei 
physical  description  and  a sketch  of  each  model  is 
given.  The  numerical  values  of  the  masses,  dampers, 
springs,  loads  and  initial  conditions  are  presented 
so  that  the  problem  can  be  modeled  by  the  Interested 
reader. 

FOUR-WHEEL  DOLLY 

The  first  model  is  a four-wheel  dolly  that  has  a 
short  transient  period  followed  by  damped  free  vibra- 
tion. The  model  is  shown  in  Figure  A-l  and  the  prop- 
erties data  are  given  in  Table  A-l.  Masses  1-4  repre- 
sent the  four  wheels,  mass  5 the  translation  inertia 
and  masses  6 and  7 the  rotation  inertias  of  the  dolly. 
Springs  and  dampers  1-4  represent  a linear  suspension 
system  and  sprlnga  5-8  represent  the  nonlinear  spring 
characteristics  of  a wheel  in  contact  with  the  ground, 
l.e.,  linear  in  compression  and  zero. force  in  tension. 

A constant  force  of  -5151. 04N  (-1158  lbs)  is  applied 
to  mass  5 to  represent  the  deed  load;  along  with  ini- 
tial dlaplacements  of  *1  * *2  “ x3  “ *4  * -0.007353m 
(-0.2895  in.)  and  xj  • -0.022060m  (-0,8685  in.)  to 
place  the  dolly  in  equilibrium.  The  excitation  is 
provided  by  an  Isosceles  triangle  pulse  of  2224. UN 
(500  lbs)  peak  at  0.025s  applied  to  mass  1,  represent- 
ing a severe  bump  encountered  by  wheel  1,  This  tran- 
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slent  load  Is  severe  enough  to  first  lift  wheel  1 off 
the  ground  and  then  wheel  4 on  the  rebound.  So  the 
notion  of  the  dolly  is  a rotation  about  the  axis  formed 
by  wheels  2 and  3 along  with  some  translation.  After 
the  transient  and  nonlinear  rebound  phase  the  dolly 
rocks  back  and  forth  with  the  motion  dampening  out. 
Hence  during  this  problem  there  Is  a transient  phase 
of  higher  frequency  response  followed  by  a long  free 
vibration  phase  of  lower  frequency. 

DROP  TEST 

The  second  model  Is  a drop  test  package  that  Is 
dropped  onto  a nonlinear  spring,  representing  the 
ground.  The  model  Is  shown  In  Figure  A-2  and  the 
properties  data  are  given  in  Table  A-2.  Hass  7 Is 
seen  to  represent  a large  rigid  portion  of  a structure 
with  Internal  packages  (components)  represented  by 
masses  1-6,  and  dampers  and  springs  2-7.  Spring  1, 
tabulated  In  Table  A-2,  is  a nonlinear  spring  giving 
an  approximate  system  frequency  of  10  Hz  that  repre- 
sents the  ground  or  some  surface  upon  which  the  pack- 
age Is  dropped.  The  first  two  oscillators,  mj  and  m2, 
are  1.0  Hz  oscillators  (undamped)  and  with  the  first 
oscillator  damped  at  102  of  critical  and  the  second 
at  1102  of  critical.  The  next  pair  of  oscillators, 
m)  and  m4,  are  10  Hz  oscillators  with  the  same  damp- 
ing characteristics  and  the  last  pair  of  oscillators, 
mj  and  mg,  are  100  Hz  oscillators  with  the  same  damp- 
ing characteristics.  The  Initial  displacement  of  all 
the  masses  (degrees  of  freedom)  Is  0. 0254m  (1.0  In.). 

A force  of  -1717. 01N  (-386.0  lbs)  is  applied  to  masses 
1-6  and  a force  of  -171701. 37N  (-38600.0  lbs)  is  ap- 
plied to  mass  7;  this  represents  the  gravitational 
force.  With  these  conditions  the  response  Is  composed 
of  periods  of  free  fall,  a shock  excitation,  rebound, 
free  fall,  etc. 


CANTILEVER  BEAM 

The  third  model  Is  a cantilever  beam  with  both 
linear  and  nonlinear  stiffness  characteristics  sub- 
jected to  a lateral  tip  load  (step  time  function). 

The  model  Is  shown  In  Figure  A-3  and  the  properties 
data  are  given  In  Table  A-3.  The  cantilever  beam  Is 
discretized  with  14  equal  mass  elements  with  rotation- 
al and  translational  degrees  of  freedom.  The  stiff- 
ness is  characterized  by  bending  and  shear  springs; 
no  damping  Is  considered.  The  discrete  elements  are 
based  on  the  properties  for  a continuum  beam  with  the 
dimensions  and  properties  as  shown  In  Figure  A-3. 

The  nonlinear  model  considers  nonlinear  bending  stiff- 
ness only;  the  shear  stiffness  remains  linear.  The 
connection  to  physical  reality  is  tenuous,  but  It  can 
be  thought  of  as  a two-strip  cantilever  beam  of  the 
same  total  thickness  that  has  been  riveted  together. 
The  moment-rotation  relationship  remains  linear  for 
a small  rotation  and  then  Is  nonlinear.  Indicating 
some  slipping  of  the  riveted  Joints.  The  purpose  of 
the  nonlinear  model  is  to  compare  the  behavior  of  the 
time  step  selection  logic  to  that  for  the  linear  model 
and  to  determine  the  computational  overhead  for  a non- 
linear model.  The  excitation  Is  a step  (at  zero  time) 
load  of  889. 64N  (200  lbs)  applied  at  degree-of- freedom 
2 In  the  positive  direction. 

AXIAL  BAR 

The  fourth  model  is  an  axial  bar  subjected  to 
wave  propagation.  The  model  and  properties  are  given 
in  Figure  A-4.  Note  the  properties  are  chosen  to 
give  a maximum  frequency  of  200  rad/sec,  hence  a maxi- 
mum time  step  of  0.01  sec  for  the  central  difference 
method.  The  Initial  condition  Is  a displacement  of 
0. 0254m  (1.0  in.)  at  degree-of-freedom  11,  i.e.,  at 
the  center  of  the  bar. 
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Fig.  A-2  Drop  test  model 
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Fig.  A-l  Four-wheel-dolly  model 


Id  '■qua l naan  eltNaeatm  of  0 Oil?  ■ <0  1 ml  leagth 
Traaalattnaal  aad  rotation. I inertia 


Fig.  A-3  Cantilever  bens  model 
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Table  A-l  Structural  dynamic  properties  of  the  four-wheel-dolly  model 


Mass 

Dampers 

Springs 

Number 

(kg) 

((lb-sec2/ln)) 

(N-S/m) 

((lb-sec2/ln)) 

(M/m) 

((lb/in)) 

1 

8.7563 

(0.05) 

700.51 

(4.0) 

87563.43 

( 500.0) 

2 

8.7563 

(0.05) 

700.51 

(4.0) 

87563.43 

( 500.0) 

3 

8.7563 

(0.05) 

700.51 

(4.0) 

87563.43 

( 500.0) 

A 

8.7563 

(0.05) 

700.51 

(4.0) 

87563.43 

( 500.0) 

S 

525.3804 

(3.00) 

175126.85 

(1000.0)* 

6 

10507.6080 

(60.00) 

175126.85 

(1000.0)* 

7 

10507.6080 

(60.00) 

175126.85 

(1000.0)* 

8 

175126.85 

(1000.0)* 

* Compression  only,  no  spring  In  tension 


Table  A-2  Structural  dynamic  properties  of  the  drop  test  model 


Mass 

Dampers 

Springs 

Number 

(kg) 

((lb- 

sec^/ln)) 

(N-S/m) 

((lb-sec/ln)) 

(N/m) 

((lb/in)) 

1 

175.1268 

( 

1.0) 

220.07 

(1.25664) 

6913.73 

(39.4784) 

2 

175.1268 

( 

1.0) 

2420.78 

(13.8230) 

6913.73 

(39.4784) 

3 

175.1268 

( 

1.0) 

2200.71 

(12.5664) 

691372.79 

(3947.84) 

4 

175.1268 

( 

1.0) 

24207.78 

(138.230) 

691372.79 

(3947.84) 

5 

175.1268 

( 

1.0) 

22007.14 

(125.664) 

69137278.51 

(394784.0) 

6 

175.1268 

( 

1.0) 

242077.85 

(1382.30) 

69137278.51 

(394784.0) 

7 

17512.68 

(100.0) 

0.0 

(0.0) 

* 

* 

* Force-Deflection  Characteristics  of  Spring  7 


Deflection 

Force 

m 

(In) 

N 

(lb) 

-2.54 

(-100.0) 

-356027.75 

(-125000.0) 

-0.007620 

(-0.3) 

-556027.75 

(-125000.0) 

-0.005080 

(-0.2) 

-1779288.80 

(-400000.0) 

-0.002540 

(-0.1) 

-444822.20 

(-100000.0) 

-1.27xlO"7 

(-0.000005) 

-4448.22 

( -1000.0) 

0.0 

(0.0) 

0.0 

(0.0) 

2.54 

(100.0) 

0.0 

(0.0) 

( 
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Table  A- 3 Structural  dynamic  properties  of  the  cantilever  beam  model 


Number 

Mass 

Spring 

1- 25  (odd) 

2- 26  (even) 

27 

28 

4.022xl06  kg-m2  (3.56xl0~5  in-lb-sec2/rad) 

0.29947  kg  (1.71x10* 3 lb-sec2/in) 

4.022x10  ® kg-m2  (3.56x10  3 in-lb-sec2/rad) 
0.29947  kg  (1.71xlO-3  lb-sec2/in) 

6327.15  m-N  (5.6xl04  in-lb)* 

1.1208xl09  N/m  (6.4xl06  lb/in) 

12654.30  m-N  (1.12xl05  in-lb)* 

2. 2416x10 9 N/m  (1.28xl07  lb/in) 

* For  nonlinear  model  the  following  tabular  values  are  used: 


9 (rad) 

Moment** 

(m-N) 

(in -lb) 

-1.00 

-3164.21 

(-2.80056xl04) 

-0.002 

-6.9599 

(-61.6) 

-0.001 

-6.3272 

(-56.0) 

0.0 

0.0 

(0.0) 

0.001 

6.3272 

(56.0) 

0.002 

6.9599 

(61.6) 

1.00 

3164.21 

(2.80056xl04) 

**  For  spring  27  multiply  by  2.0 


ml_21  - 175.12685  kK  (1.0  lbsec2/ln) 
kl-20  " 1751268 . 5 N/m  (10000.0  lb/ In ) 
Boundary  Conditions:  Free-Free 

Fig.  A- 4 Axial  bar  model 
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